A SIMPLE FINITE DIFFERENCE METHOD FOR 
TIME-DEPENDENT, VARIABLE COEFFICIENT STOKES FLOW ON 

IRREGULAR DOMAINS 



Abstract. We present a simple and efficient variational finite difference method for simulating 
time-dependent Stokes flow in the presence of irregular free surfaces and moving solid boundaries. 
The method uses an embedded boundary approach on staggered Cartesian grids, avoiding the need 
for expensive remeshing operations, and can be applied to flows in both two and three dimensions. 
It uses fully implicit backwards Euler integration to provide stability and supports spatially varying 
density and viscosity, while requiring the solution of just a single sparse, symmetric positive-definite 
linear system per time step. By expressing the problem in a variational form, challenging irregular 
domains are supported implicitly through the use of natural boundary conditions. In practice, 
the discretization requires only centred finite difference stencils and per-cell volume fractions, and is 
straightforward to implement. The variational form further permits generalizations to coupling other 
mechanics, all the while reducing to a sparse symmetric positive definite matrix. We demonstrate 
consistent first order convergence of velocity in L 1 and L°° norms on a range of analytical test cases 
in two dimensions. Furthermore, we apply our method as part of a simple Navier-Stokes solver to 
illustrate that it can reproduce the characteristic jet buckling phenomenon of Newtonian liquids at 
moderate viscosities, in both two and three dimensions. 
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1. Introduction. The equations of Stokes flow describe the motion of fluids 
where the nonlinear advection terms present in the Navier-Stokes equations have been 
eliminated or ignored. This approximation is important in many physical scenarios 
where inertia is essentially negligible (Re -C 1), however it typically leads to the 
steady- state Stokes equations. In this paper we are instead interested in the time- 
dependent or unsteady case, motivated by its use as a building block for solving the 
Navier-Stokes equations. That is, a number of methods apply operator splitting to 
the Navier-Stokes equations in order to treat the advective terms separately from the 
time-dependent Stokes equations. 

The goal of this paper is to derive and validate a simple, efficient, and stable 
method for solving the time-dependent Stokes flow equations in the presence of irreg- 
ular free surfaces and solid boundaries, while supporting spatially varying viscosity 
and density coefficients. We solve these equations on the classic staggered Cartesian 
grid using a finite difference approach. However, in order to support irregular domains 
which arise frequently in the presence of evolving liquid interfaces and moving ob- 
jects, we must determine a proper discretization for non-grid-aligned (ie., embedded) 
boundary conditions. The free surface condition poses a particular challenge because 
it involves a delicate coupling between the boundary normal, the pressure, and the 
components of the deviatoric stress tensor. Our approach will be to take inspiration 
from the finite element method, and express the Stokes flow problem in a variational 
form relating velocity, pressure, and deviatoric stress. We will show that this yields a 
hybrid approach that is discretized with simple finite differences, but which implicitly 
enforces complex embedded boundaries through the natural boundary conditions of 
the problem. We believe this to be the first method that can provide a fully implicit 
discretization of this problem on Cartesian grids, while simultaneously yielding a sym- 
metric positive-definite linear system and correctly handling both variable coefficients 
and the difficult free surface case. 
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2. Background. While there exist finite element and finite volume methods that 
can be applied to the Stokes equations on irregular domains, these methods require 
meshes comprised of well-shaped elements that align with the physical boundaries. 
Frequent remeshing becomes necessary in settings with rapidly evolving boundaries, 
and as others have noted, such remeshing is a challenging and expensive problem in 
its own right [4]. Unstructured meshes often also incur a performance penalty as 
compared to regular grids due to the overhead costs of manipulating more complex 
data structures. We therefore pursue the simpler and more efficient embedded finite 
difference approach, allowing boundaries to cut arbitrarily through an underlying 
regular Cartesian grid. 

Viscous free surface flows have long posed challenges for finite difference methods. 
The core problem is ensuring that the region external to the liquid applies zero force 
on the surface itself through the application of appropriate boundary conditions. The 
seminal marker-and-cell (MAC) paper of Harlow and Welch [22] discussed these dif- 
ficulties, and dispensed with the free surface viscous stress conditions for simplicity. 
However, they noted that at lower Reynolds numbers the more accurate condition 
becomes necessary. Hirt and Shannon proposed a simple correction to the normal 
stress [53], which Nichols and Hirt subsequently improved to address the tangential 
stress [33] . Because an explicit discretization of viscous terms suffers from a stringent 
time step restriction of At < 0(pAx 2 / fx) that becomes problematic at high viscosi- 
ties, Pracht incorporated the preceding ideas into a fully implicit integration scheme 
to improve stability [39]. All of these schemes assume that the surface normal is 
either aligned with coordinate axes, or at a 45° angle, in order to explicitly design 
conditions for these special cases. This effectively rasterizes or voxelizes the domain 
into a "stairstep" approximation that may not realistically reflect the physical geome- 
try. This can limit the generality and accuracy of the approach, and the case-by-case 
analysis somewhat complicates the implementation. Nevertheless, Tome et al. devel- 
oped a three-dimensional extension that convincingly simulates a range of free surface 
phenomena (49] |48] , while using explicit time integration and a pressure projection 
method that splits the Stokes equations into separate pressure and viscous compo- 
nents (see [2T] for an exploration of the issues raised by this type of splitting). Oishi 
et al. later implemented a more stable implicit time integration scheme [35] within the 
same framework to enable larger time steps and improve efficiency. Aside from our 
own work, this is the only other implicit finite difference method that has addressed 
the phenomenon of viscous jet buckling with proper free surface conditions. However, 
it requires the same case-by-case analysis as its predecessors, as well as the use of 
the less robust bi-conjugate gradient method to solve a large, sparse, non-symmetric 
linear system. A few authors have also advocated a semi-implicit treatment of spa- 
tially varying viscosity, in which terms that couple different velocity components are 
treated explicitly, and non-coupled terms are treated implicitly [271 S3 EE] • 

In computer animation, staggered grid schemes have also been used extensively 
for free surface flows, again using a pressure projection approach [3 [18j [42] ■ Recently, 
Batty and Bridson showed that improper viscous free surface conditions in these 
methods tend to destroy angular momentum [3]. To address this, they proposed a fully 
implicit scheme for variable viscosity that yields symmetric positive-definite linear 
systems, and presented animations of rotating and buckling viscous liquids. While 
that method shares much with the present work, it differs in several key respects. 
First, it uses a pressure projection approach to split up the Stokes equations which 
entails weaker coupling between pressure and viscous terms. Secondly, it enforces a 
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simplified free surface boundary condition that incorrectly neglects the interaction 
between pressure and viscous stresses at the surface. Thirdly, the present approach 
yields a linear system in terms of stresses, rather than velocities. Finally, the current 
work provides numerical experiments indicating that our scheme is indeed convergent. 

There exists a broad family of methods that seek to accurately enforce boundary 
conditions of various kinds on irregular domains, while solving equations on Cartesian 
finite difference grids. We will collectively refer to this family, which includes the cur- 
rent work, as embedded boundary methods (although other authors have sometimes 
used the term immersed boundary methods [31]). As noted above, the motivation 
for embedded approaches is to avoid the substantial computational cost of generating 
and manipulating high quality, fully unstructured, boundary-conforming meshes as 
required by finite element or finite volume methods. Common examples of embedded 
boundary methods include the traditional immersed boundary methods (IBM) [38], 
ghost fluid methods (GFM) [115], immersed interface methods (IIM) matched 
interface and boundary methods (MIB) 52 , and cut-cell methods [23]. These have 
been quite effective for a number of problems and some can achieve higher order 
accuracy. Nevertheless, we are not aware of any symmetric fully implicit embedded 
boundary method with a comparably simple implementation that can address the 
variable coefficient unsteady Stokes equations with free surfaces. As an example, re- 
cent work that primarily emphasizes the use of the ghost-fluid approach for simulating 
turbulent atomization [13] . Desjardins et al. revert to a diffuse interface method for 
viscous terms, stating that this choice was motivated by the complexity of viscous 
GFM discretizations [55] and the difficulty of achieving an implicit formulation. 

Among existing embedded boundary methods, the most closely related to our 
approach are staggered grid methods that use standard centred-difference stencils, 
scaled by simple diagonal weighting matrices to support irregular boundaries. These 
methods yield symmetric positive-definite linear systems, but have primarily been 
applied to Poisson and diffusion problems to date. First, a simple ghost fluid method 
has been used to enforce Dirichlet boundary conditions [8j [2Ql [TTl [33] for pressure 
projection methods and for implicit integration of spatially constant viscosity with 
solid boundaries. Secondly, a finite volume-like technique has been proposed for han- 
dling Neumann (and Robin) boundary conditions in the context of pressure projection 
methods for solid-fluid interaction [40j [45j [33l [37] . Work by Batty et al. on solid-fluid 
coupling and viscous flows E] also falls into this category, although they are de- 
rived from variational principles rather than ghost fluid or finite volume concepts. 
The current work directly extends and unifies these last two methods, while provid- 
ing numerical experiments that validate the approach. 

Robinson-Mosher et al. have developed a related symmetric positive-definite for- 
mulation of time-dependent Stokes flow, focusing on solid-fluid coupling [43]. How- 
ever, their viscosity discretization is limited to voxelized (axis-aligned) solid geometry 
and constant viscosity, and does not address the free surface boundary condition 
we consider. Moreover, our derivation identifies their algebraic transformation as 
a change of variables from a velocity-pressure formulation to a pressure-stress for- 
mulation, lending physical insight to the method, and recovers a more fundamental 
variational form of the mechanics that offers improved specialized solvers. 

There arc of course a host of methods for treating the more common symmet- 
ric indefinite form of the Stokes equations for steady and unsteady problems [5]. 
Among these are several methods carefully designed to be optimal for the Stokes 
problem, in the sense that their convergence rates are independent of the size of the 
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discrete simulation mesh, with "appropriate" choices of preconditioners or smoothers 
[Tl] . Nonetheless, such methods face additional challenges in addressing boundary 
conditions and variable coefficients, particularly for the case of irregular embedded 
boundaries. For example, Wan et al. applied substantial modifications to a geometric 
multigrid method to efficiently handle embedded interfaces for the simpler Poisson 
equation [51] . Similarly, many fast solvers for the Stokes problem rely on a simplifi- 
cation (discussed in section 13. fj) that decouples the viscous contributions into three 
independent Poisson-type problems, enabling the use of existing fast Poisson solvers 
as sub-components. This simplification does not readily apply in the free surface 
or variable viscosity settings. We therefore prefer to pursue a symmetric positive- 
definite formulation for which effective black-box iterative solvers and preconditioners 
are more readily available. The development of special-purpose solvers and precondi- 
tioners for the positive-definite formulation is a potentially exciting research direction, 
which we defer to future work. 

3. Time-Dependent Stokes Flow. The governing equations for time-dependent 
Stokes flow are: 

pu t = -Vp + V -t (3.f) 
V-u = (3.2) 
„ (Vu+{Vu) T \ 
T = 2M [ 2 J (3 - 3) 

where p is the fluid density, u is the fluid velocity, p is the pressure, r is the symmetric 
deviatoric stress tensor, p is the dynamic viscosity coefficient, and the subscript t indi- 
cates a time derivative. We allow both p and p to vary spatially. At solid boundaries, 
the no-slip condition applies, which dictates that the fluid velocities match those of 
the solid boundary, ubc'- 

u = ubc (3.4) 

At a free surface, the fluid is subject to the constraint that the traction T applied at 
the surface is zero: 

f={-pI + T)n = 6 (3.5) 

In the above, ft is taken to be the outward normal of the free surface and / represents 
the identity matrix. 

3.1. A Note on the Simplified Stokes Equations. In situations where the 
viscosity is spatially constant, some manipulations are often carried out to reduce the 
contribution of viscosity to a simpler form. Specifically: 

V • r = V ■ (p(Vu + (Vu) T )) (3.6) 

= pV-Vu + pV- (V5) T (3.7) 

= pV - Vu + ^V(V ■ u) (3.8) 

= pN ■ Vu (3.9) 

where a simple vector calculus identity and the divergence free property of u have 
been used to eliminate the second term on the right hand side. This reduces the 
contribution of viscosity to a simple Laplace operator applied to each component of 
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velocity independently. However, the natural boundary conditions of this modified 
problem are quite different from those of the original problem, as discussed by Limache 
et al. [25J if care is not taken this can lead to non-physical solutions, particularly 
for free surface flows. We therefore prefer the fully general form of the viscous terms, 
even for constant viscosity [3j. 

4. Variational Formulations for Time-Dependent Stokes Flow. A back- 
wards Euler time discretization of the governing equations yields: 

-L(u-iT) = V-T-Vp (4.1) 

V-u = (4.2) 
r = /i(W+ (Vu) T ) (4.3) 

In these expressions u is the final velocity field, u* is the input (possibly divergent) 
velocity field, and At is the size of the time step. 

An equivalent variational formulation is the following: 

Pi,.-, .-*..2 a,„vt , a,_ J^7u+{Vu) T \ A* 2 



max nun - u*\\ 2 - AtpV ■ u + Atr : ^ ^ >- j - — ||r|| 2 F (4.4) 

The domain of integration is the liquid (non-air) region, VLl. Calculus of variations 
can be used to show that the optimality conditions for this problem yield precisely 
the equations of the original PDE problem above, with the zero traction free surface 
boundary condition (|3.5|) enforced as a natural boundary condition. 

Similarly, the following variational formulation enforces the Stokes equations in- 
side the domain, but with static solid boundary conditions (u = 0). 

maxmin /// -\\u - u*\\ 2 + Atu ■ (Vp - V ■ r) - — ||r||| (4.5) 
p> t « JJJqf 2 4/i 

Here the integration is performed over the fluid (non-solid) region, f2 p . Note that these 
two variational formulations differ essentially by an integration by parts operation on 
the middle terms. 

5. Discretization. With a particular variational formulation in hand, we pro- 
ceed to directly discretize the required integrals in a manner similar to that proposed 
by Batty et al. [HO]. This results in a discrete optimization problem in which the 
boundary conditions are enforced naturally and implicitly. This is in contrast with the 
more common approach of first discretizing the PDE form itself, which necessitates 
the explicit handling of the potentially difficult boundary conditions outlined earlier. 

We discretize the derivatives using centred finite differences on the classic stag- 
gered (MAC) grid, with pressures at cell centres, and velocity components on cell 
faces [22] . To support viscosity we must also place the components of the deviatoric 
stress tensor on our grid. The most natural way to do this is to locate diagonal com- 
ponents of the stress tensor (t xx , T yy , t zz ) at cell centres, and off-diagonal components 
T yz ) on cell edges (nodes in two dimensions). This arrangement is illustrated 
in 2D and 3D in Figures 15.11 and 15.21 respectively. Straightforward centred differenc- 
ing can then be used to compute the required derivatives in the correct locations. 
This approach was first proposed by Darwish et al. |12j for two dimensions, and later 
extended to three dimensions by Mompean and Deville [32] . While it has traditionally 
been applied to non-Newtonian fluids in which the constitutive equations are more 



G 



FINITE DIFFERENCE STOKES FLOW ON IRREGULAR DOMAINS 



□ — I — H — I — H — I — H — I — m 



[] — I — E3 — I — E3 — I — E3 — I — E] 



O 



O 



O 



O 



O 



O 



[] — I — B — I — B 



O 



O 



O 



O 

— I— 

O 



[] — I — E3 — I — E3 — I — E3 — I — E] 



O 



O 



O 

ffl — I — Bp 



O 



O 



□ — I — B — I — B — I — B — I — B 

Fig. 5.1. The standard staggered pressure-velocity grid layout in 2D, with stress components 
added. Circles indicate the locations of pressure and diagonal stress components (t xx )■ Dashes 
across cell faces indicate horizontal and vertical velocity components. Small squares at cell corners 
indicate the locations of off-diagonal stress components (r xy ). 



complex, we find its simplicity and elegance useful for the purely Newtonian flows we 
consider. 

Note that because the deviatoric stress r does not measure compression, we are 
assured that Tr(r) = which in 2D implies that t xx = —t vv . We can therefore simplify 
the necessary computations by solving just for t xx rather than both quantities; in 
the final linear system, this will yield equations of the form t xx = ft — , 
thereby retaining symmetry. Similar transformations apply in 3D to eliminate t zz — 
-(t^ + T yy ), yielding t xx + \r yy = fi (ff — §7) and r yy + \t xx = fj, (§H - |^J. 
Naturally, the stress tensor will also be symmetric, so that r xy = r yx , t xz = t zx , and 
T yz = r zy , which further reduces the number of variables to be computed. 

We approximate the integral comprising the variational forms by scaling each term 
by the fractional volume of material in a cell-sized control volume surrounding the 
appropriate sample point, and sum over all cells. For example, terms that lie on faces 
are scaled by the volume of fluid in a square (or cubic) control volume surrounding the 
face centre. The necessary two-dimensional control volumes are illustrated in Figure 
15.31 (While we could apply a higher order approximation of the integrals, it is this 
simple piecewise constant approximation that ensures we retain the same stencils as 
a standard grid-aligned finite difference discretization.) For free surface boundary 
conditions, we need to estimate volume fractions interior to the liquid (ie. not air), 
indicated by weights Wl. Later, we will also need the corresponding air fractions, 
Wa — I — Wl (where / the identity). Likewise, for solid boundary conditions we 
estimate the volume fractions Wf of a cell that is inside the fluid (ie. not solid) , and 
its complementary solid fraction, Ws = I—Wf. These volume fractions are illustrated 
in Figure [5T4l and can, for example, be estimated from a level set representation using 
the method of Min and Gibou [30]. We will see that this simple piecewise constant 
approximation of the integrals ensures that we retain simple centred finite difference 
stencils, while volume fraction weighting handles boundary conditions. 

In equation (|4.4[) . we are integrating over the liquid region f^, so we must use 
volume fractions with subscript L. The first term of the equation consists of velocity 
data that lies on faces, so we estimate the integral using fractional volumes associated 
to each velocity face sample. We indicate this using superscripts on the weights that 
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Fig. 5.2. The standard staggered pressure-velocity grid layout in 3D, with stress components 
added. The black circle indicates the location of pressure and diagonal stress components (r X x,T yy ). 
The colored squares on cell edges indicate the locations of off-diagonal stress components (r yz is red, 
t~xz ^ green, T xy is blue). Dashes across cell faces indicate velocity components (u is red, v is green, 
w is blue). 
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FlG. 5.3. Control volumes around each sample location in 2D. From left to right: pressure 
and diagonal stress control volume, off-diagonal stress control volume, horizontal velocity control 
volume, vertical velocity control volume. 



indicate the type of control volume being considered; in this case the weight used 
is Wf^. Similarly, the second term consists of pressures and divergences that are 
conceptually located at cell centres, so we use cell-centred weights, W^. Finally, the 
last two terms are associated with stresses, which lie on cell edges and centres, so we 
use the fractions W£. 

Applying this framework to (|4.4I) leads to the following discrete Stokes problem 
with free surface boundaries: 

maxmin -{u - u*) T PWf){u - u*) + Atp T W?G T u + Atr T W£Du - —t t M~ x W t l t 

(5-1) 

In this expression P is a diagonal matrix of densities per velocity sample and M 
is a diagonal matrix of viscosity coefficients per stress sample. The derivatives are 
discretized with staggered grid finite difference operators: G is the usual discrete 
gradient, and D is the discrete deformation rate applied to velocity (ie. Du w 
|(V«+(Vti) T )). Note that the negative transposes of G and D are the corresponding 
discrete vector and tensor divergence operators, respectively. The various W terms 
are diagonal matrices consisting of the volume fractions described above. 

Since this discrete optimization problem is a quadratic, the optimality conditions 
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Fig. 5.4. An illustration of volume fraction regions for solid and free surface weighting. Left: 
A geometric scenario in which a solid (dark gray) meets a body of liquid (blue) in the presence of 
air (white). Note the air-liquid interface extrapolated into the solid. Middle: The volume fraction 
region for a liquid (ie. non-air) weight, Wr,, shown in gray, and its complementing air fraction Wa 
shown in white; the presence of the solid is ignored. Right: The volume fraction region for a fluid 
(ie. non-solid) weight, Wf, shown in gray, and its complementing solid fraction, Ws, shown in 
white; in this case the position of the liquid-air interface is ignored. 



yield the following symmetric indefinite linear system: 

H^l DT Wi GWl\ ( u\ ( i- t PWtu* 
W[D -\M~ l Wl t = | (5.2) 

W P L G T J \ p J \ 



The solid boundary Stokes problem (|4.5j) can be discretized in much the same 
manner, except that integrals are computed over the fluid (non-solid) region using 
volume fractions Wf- The discrete form is: 

maxmin i(w - u*) T PW^{u - u*) + Atu T W^{Gp + D t t) ~ ^t t 'M -1 W£t (5.3) 
p,t u 2 4 

The linear system for the solid wall Stokes problem is: 

^PW% W^D T W^G \ ( u\ ( ±PW£u* \ 
DW£ -\M- X W^ r = (5.4) 

G T Wp / \ / 

Because these two systems have nearly identical forms, and differ only by the 
weighting matrices W, we can straightforwardly combine them to handle both free 
surfaces and solid boundaries in the same problem: 

■^.PWpWl WfD T Wl W^GWl \ I u \ / ±PW%WZu* 
W[DWZ -\M- x W T F Wl II t 1 = 1 
Wf j G T Wp I \ P J \ 

(5-5) 

By combining the two formulations only at the discrete level, we are able to exploit 
the natural boundary conditions to handle both boundary types together, despite the 
fact that in each of the two continuous variational formulations only one of the two 
boundaries can be considered natural. 

For this Stokes system, eliminating stress by applying the Schur complement to 
the centre block yields the symmetric indefinite velocity-pressure system most com- 
monly associated with the Stokes problem. However, as outlined earlier we prefer 
to solve symmetric positive-definite systems, as they are generally more amenable to 
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common black box solvers such as preconditioned conjugate gradient methods, domain 
decomposition, etc. Conveniently, the upper-left block of the full system is a diagonal 
matrix which can be trivially inverted. We therefore perform a Schur complement 
on this block to eliminate velocity and arrive at a sparse symmetric positive-definite 
system for pressure and stress. This is exactly analogous to classic pressure projec- 
tion methods for incompressibility, in which eliminating velocity yields a symmetric 
positive-dehnite Poisson equation for pressure alone. Of course, it should be em- 
phasized that this technique can only be applied in the time-dependent Stokes case; 
otherwise the upper-left block is simply zero making a Schur complement impossible. 
The symmetric positive-definite form is: 

f An A 12 \ ( t\ _( W[DW£u* 
\AJ 2 A 22 ){p )-{ WiG T W£u* 

where the blocks of the matrix are 

A n = ^M^WlWp + AtW[DP- 1 W^ 1 WpD T Wl 

Au = MWlDP- l WZ~ 1 W%GWl 
A 22 = AtW p L G J ' P- 1 W£- 1 W$GWZ 

Our results for the Stokes problem consistently indicate hrst order convergence 
for all variables in L , and hrst order convergence for velocity in L°° . Stress fails 
to converge in L°° due to noisy errors along the boundaries. Given that stress is 
computed as the gradient of velocity, it is not entirely surprising that it loses one 
order of accuracy in L°°. Fortunately, for most practical applications the behaviour 
of velocities is of greater importance. 

6. Non-Homogeneous Boundary Conditions. The preceding formulation 
for the Stokes problem exploits natural homogeneous boundary conditions to simplify 
handling of irregular domains on Cartesian grids. However, many practical situa- 
tions will call for non-homogeneous boundary conditions where the boundary values 
are non-zero. For example, moving solid boundaries will require non-zero boundary 
velocities to be enforced, as will inflow and outflow boundaries. Similarly, non-zero 
pressure boundary conditions have been used to support surface tension (eg. [T7]). 
To incorporate non-homogeneous boundary values into our framework in a consistent 
manner, we introduce additional terms to account for the work done by the boundary 
itself. 

6.1. Prescribed Traction Boundaries. To add a prescribed traction bound- 
ary, we need to account for the work done by the traction at the surface. We do this 
by adding the following boundary term to (I4.4[) : 

// Atu ■ (pbcI ~ TBc)n (6.1) 
J J an A 

where pbc an d tbc are the prescribed pressure and deviatoric stress, respectively, 
and n in this case is the outward normal with respect to the air (non-liquid) region, 
Qa- (Of course, only the resulting normal component, ie. surface traction, will 
be enforced.) This surface integral can be converted to a volume integral through 
integration by parts: 

At JJJ PBcV ' U ~ TBC : ( J U+ ^ U) ^j+u- (VpBC-y -tbc) (6.2) 
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This ensures that all the relevant quantities are at grid locations that are consistent 
with the functionals considered previously. This also avoids the need to discretize 
surface integrals which would add further complications and potentially sacrifice the 
convenience of centered difference stencils. 

Using Wa terms to indicate the air fraction of a particular control volume, the 
discretized form is: 

At {- P T BC W p A G T u - tI c W t a Du + u T Wl{G PBC + D t t bc )) (6.3) 

The new terms modify the right hand side of the linear system (I5.2j) . to become: 

£PW£u* + GWIpbc ~ WIGpbc + D t W t a t bc - W u a D t t bc 

j (6.4) 



Although the air region may extend far from the actual liquid surface, we only need 
to apply modifications to the right hand side for rows of the system in which the 
matrix has non-zero entries, indicating that there is liquid present. This means that 
in practice only the interface between the air and liquid regions plays a role. 

6.2. Prescribed Velocity Boundaries. In the common case of moving solid 
boundaries with prescribed velocities ubc > we account for the work done by the solid 
on the fluid by adding the following term to (14.51) : 

Atu B c ■ {pi — r)n (6-5) 

dn s 

where n is the outward normal to the solid region, £1$. In volume integral form we 
have: 

*tjjf a P V ■ ubc r : ( VgflC ±iYggg£ ) + Ubc ■ (Vp - V • t) (6.6) 

Labelling solid fractions Ws and discretizing, we arrive at the following term: 

At {-p T W p s G T u BC - t t W t s Du BC + u T BC Ws{Gp + D t t)) (6.7) 

This results in a modification to the right hand side of the linear system (I5.4[) . to 
become: 

-tt PW >* 
-DW^ubc + W T s Du BC 
-G T W%u BC + W p s G T u BC 

These right hand side modifications are also only applied to rows in which the matrix 
has valid non-zero entries. 

7. Generalization and Further Factorization. At this point, we highlight 
the abstract form of the problem studied above, which helps see the discrete na- 
ture of the transformation we use and illustrates how to generalize our technique to 
other problems such as two-way fluid-solid interaction. This is a special case of a 
weighted linear least-squares problem with block diagonal regularization subject to 
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linear equality constraints, which we state using bold elements to avoid confusion with 
prior notation: 

min i(u-u*) T M(u-u*) + i(b- Au) T W(b - Au). (7.1) 
u 

Cu + d = 



For the Stokes discretization above, the velocity vector u here corresponds to the 
vector u of all fluid velocities and M, the block diagonal regularization matrix which 
we term the mass matrix for the abstract problem, corresponds to the (simply di- 
agonal) product PWpW^. The first term in the objective is thus a kinetic energy 
norm of the difference between the old (or predicted) velocity u* and the new veloc- 
ity u which we are finding. The constraint matrix C and constraint vector d cor- 
respond to the negative weighted discrete divergence operator AtW\G T Wp and the 
zero vector (for the divergence of the Stokes velocity in the interior) modified as neces- 
sary by non-homogeneous boundary conditions. Finally, in the weighted least-squares 
term the matrix A corresponds to DWp (i.e. the weighted measure of deformation 
rate), the vector b is zero in the interior of the Stokes problem but modified by non- 
homogenous solid boundaries as needed, and the weighting matrix W corresponds to 
2AtW£MW F - 1 (i.e. the weighted, time-scaled viscous coefficients). 
Rearranged into this format, our Stokes discretization looks like: 

min ^(u-u^PW^WKu-u^ + AtiDW^ufWlMW^iDW^u) 

u 

AtW p L G T W^u = 

(7.2) 

The solution is a balance between staying close to the old velocity (weighted by 
density) and minimizing the deformation rate (weighted by viscosity and time step) 
with cell fractions accounting for irregular geometry, subject to incompressibility. This 
makes evident the connection to Helmholtz's minimum dissipation theorem for Stokes 
flow pQ. 

As before, we can write down the KKT optimality conditions, introducing a 
Lagrange multiplier A for the linear constraint: 

M + A T WA C T \ / u \ _ / Mu* - A T Wb \ , , 

c o )\ \ )-\ -d )■ [7 - d > 

Up to a scale factor, the new variable A corresponds to pressure in the Stokes problem. 

This is of course a symmetric indefinite matrix, but our variable change to pressure 
and viscous stress generalizes here as well to arrive at a sparse SPD matrix. First 
introduce the weighted least-squares residual 

r = W(b-Au), (7.4) 

which in the Stokes problem is the weighted and time-scaled viscous stress. Augment- 
ing system (|7.3I) with equation (|T.4[) multiplied by W 1 we arrive at: 



(7.5) 
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Eliminating velocity with the first equation u = u* + M 1 A T r — M 1 C r A leaves us 
with the SPD Schur complement problem for the residual and the Lagrange multiplier: 



W^+AM^A 1 -AM _1 C T \ /r\ /b-Au 
-CM^A 7 CM-'C T { A ~[ d + Cu 



(7.6) 



A key advantage to this system is that, under the assumption that the mass matrix M 
and the weight matrix W are block diagonal or more specifically have sparse inverses, 
and that A and C don't have overly dense columns, this is a sparse matrix itself. (If 
W isn't of this form, using f = b Au can restore sparsity.) 

Being sparse and SPD, many black-box linear solvers are available which may not 
apply to the original indefinite system. In our studies, we used PCG with a parallelized 
algebraic multiplicative Schwartz overlapping domain-decomposition preconditioner, 
for example. 

We can go further, however, reducing the requirement on A and C. In particular, 
the matrix in equation (|7.6p can be factored conveniently as: 



W-'+AM-'A 7 -AM"'C T 
-CM~'A T CM^C T 



-A I \ / M- 1 \ / -A T C T 

C O/l w- 1 Hi 



(7.7) 



These factors are obviously as sparse as the original components. Note also that the 
right-hand side of equation (|7.6p can be written as 



b - Au* \ _ f -A I \ f M- 1 \ ( Mu* \ / 

d + Cu* ) ~ [ C oil W- 1 [ Wb j + I d 



(7.6 



In the case of Stokes with homogeneous boundary conditions, d = 0, and it is then 
clear that equation (|7.6j) gives the normal equations for a sparse, weighted, and un- 
constrained least-squares problem: 



mm 



Mu* \ _ / -A T C T \ ( r 

wb [ i o Ma 



(7.9) 

[M-iW- 1 ! 



Specialized solvers such as LSQR [36) may then offer a significant advantage. 

The primary advantage of this generalization, and factorization, lies however in 
the ease of extending the formulation to more general dynamics. For example, Batty et 
al. 2 coupled rigid body dynamics with inviscid incompressible flow using a pressure- 
only subset of equation (|7.6[) expressed as an unconstrained least-squares problem. 
However, for rigid bodies overlapping many fluid grid cells the system suffered a large 
dense block corresponding to those cells, leading Robinson-Mosher et al. to pursue an 
indefinite form [44) . In the factored form, however, the rigid body merely corresponds 
to six fairly dense rows in A and C, which can be handled much more efficiently in 
both storage and multiplication with the factored form. Elastic and/or constrained 
solid dynamics can be phrased in the same minimization form as equation (|7.ip , 
e.g. English and Bridson's isometrically deforming membranes [15] , further opening a 
clear route to efficiently solvable general solid-fluid coupling: the additional degrees 
of freedom of the solid arc appended to u, the solid mass matrix is appended to M, 
additional solid constraints are appended to C and d, and elastic potential energy 
terms are expressed as a possibly non-linear least-squares terms appended to A and B 
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Fig. 8.1. A case in which a null space arises for pressure in a free surface problem. The 
top right cell contains some liquid, and therefore will have a positive volume weight associated with 
the divergence constraint on the cell. However, two of its associated velocity face control volumes 
(indicated by dashed blue squares) contain no liquid and will have zero volume weights. We identify 
this pressure sample as invalid, and replace it with the value of the boundary condition (eg. p = 0), 
thereby eliminating the null space. 

(as is quite natural for hyperelastic finite element models, and points out the utility 
of a Hessian-free Gauss-Newton iteration for solving them). Finally, we note that 
while Robinson- Mosher et al. discuss some of the above manipulations |33], they do 
so for a less general linear algebra problem, and do not discuss the optimization form 
for either the initial constrained problem nor the final unconstrained least-squares 
problem. 

8. Null Space Elimination. An issue that arises with our approach is the 
presence of null spaces due to overlapping volume weights assigned to different terms 
of the discrete variational problem. For example, consider the divergence operator 
for a cell near a free surface boundary; there are volume weights associated with 
both velocities (face centres) and pressures (cell centres). Cases frequently arise in 
the discretized system in which a pressure with a non-zero volume weight enforces a 
divergence constraint on at least one velocity face with zero associated volume. This 
velocity sample will appear in no other equations because of its zero volume weight, 
and therefore it may take on an arbitrary value as long as it satisfies the constraint. 
An example of such a null space scenario is shown in Figure 18.11 

In the final result, only samples with a positive associated volume weight are 
considered, so the physical solution is not adversely affected. However, such large 
spurious null spaces can pose problems when applying standard solvers for sparse 
linear systems; we would therefore like to eliminate them. 

To do so, we identify each variable that enforces a relationship on a sample 
with zero volume weights. For example, in the free surface pressure case a non- 
zero weighted pressure is tagged as invalid if it enforces the divergence-free condition 
on one or more velocity faces with zero weights. Likewise, in the solid boundary 
case if a non-zero weighted velocity sample borders a zero-weighted pressure sample, 
that velocity sample is tagged as invalid. This process can likewise be extended to 
the viscous terms. Once these invalid variables have been identified, they can be 
straightforwardly eliminated from the linear system, and replaced with the value of 
the boundary condition, resulting in a modification to the right-hand-side. This re- 
duced set of equations retains symmetry and has the same solution, but no longer 
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Table 9.1 
Solver Behaviour 





PCG Iterations 


PCG Iterations 


Grid 


Free Surface IpTl) 


Solid Wall 


16* 


1 


1 


32 2 


1 


1 


64 2 


6 


4 


128 2 


12 


14 


256 2 


27 


42 


512 2 


54 


91 


1024 2 


110 


195 



suffers from large null spaces. 

There can be one additional null space when the fluid domain is completely en- 
closed by solid walls. Because only the gradient of pressure affects the final velocities, 
pressure solutions that differ by a constant are effectively equivalent. This rank 1 null 
space doesn't pose substantial problems, so we have not eliminated it; if necessary, it 
could be removed by arbitrarily fixing one pressure value in the domain. 

9. Convergence Studies. We have verified that the method computes the exact 
solution for linear problems even for irregular domains. This includes the case of 
hydrostatic fluid with solid (and possibly free surface) boundaries, as well as rigid 
translations and rotations of liquid bodies with free surface boundaries. We provide a 
range of examples below to illustrate the convergence orders achieved by our methods 
in more difficult scenarios. All of the examples make use of curved boundaries which 
do not align with the underlying Cartesian grid, and we consider a single time step 
of the time-dependent problem in question. The examples test our methods in the 
presence of free surfaces and both static and moving boundaries. We compute the 
L°° and L 1 errors, where our discrete L 1 norm is computed as = J2i \ u i\h> f° r 

a uniform grid spacing h = Ax in d spatial dimensions. 

For each case, we transformed the linear system to the sparse symmetric positive- 
definite form as described above, and solved it with the conjugate gradient method 
preconditioned with overlapping multiplicative Schwarz domain decomposition. The 
preconditioner was determined purely algebraically, from a simple graph partition of 
the sparse matrix; we expect the positive-definite form would allow the use of many 
other black-box solvers as well. Though we do not consider this iterative solver a core 
contribution of the paper, Table l9Tl presents some representative iteration counts to 
illustrate its scaling. 

9.1. Stokes Flow with Free Surface Boundaries (2D). Our free surface 
Stokes test case is a fluid disk of radius r = 0.75 centred at the origin, with density 
p = 1 and viscosity [i = 0.1, computed over a timestep At = 1 . For simplicity 
of presentation, we describe the final velocity field in terms of a streamfunction, tp, 
where the velocity field can be derived as u final = V x tp. This also guarantees that 
the velocity field is divergence free. The streamfunction is: 

128 

4> = — r 4 cos(26») cos(v / 3 In r) (15 - 30r + 16r 2 ) (9.1) 
81 

This is a non-trivial velocity field designed to fulfill the free surface zero traction 
condition at r = 1, smoothly blended into a zero velocity at the origin (r = 0). The 
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Table 9.2 

Convergence of Stokes with free surface ( 2D ) 



Grid 


\\p- p"lloo 

1 1 if f || oo 


Order 


\\r F 111 


Order 


16 2 


1.2428E-001 




1.4886E-001 




32 2 


1.9439E-001 


-0.65 


5.3215E-002 


1.48 


64 2 


1.6330E-001 


0.25 


2.2135E-002 


1.27 


128 2 


1.3291E-001 


0.30 


8.1145E-003 


1.45 


256 2 


1.7242E-001 


-0.38 


5.3320E-003 


0.61 


512 2 


1.7261E-001 


-0.00 


2.1992E-003 


1.28 


1024 2 


1.6894E-001 


0.03 


1.1239E-003 


0.97 


Grid 


llr - T h II 

1 1 xx XX 1 1 00 


Order 


| r TT — riLlh 

II XX XX II 1 


Order 


16 2 


2.7154E-001 




1.7116E-001 




32 2 


1.4558E-001 


0.90 


3.8890E-002 


2.14 


64 2 


7.8877E-002 


0.88 


1.4642E-002 


1.41 


128 2 


5.5237E-002 


0.51 


5.6140E-003 


1.38 


256 2 


7.0687E-002 


-0.36 


3.2991E-003 


0.77 


512 2 


7.1097E-002 


-0.01 


1.3680E-003 


1.27 


1024 2 


7.4305E-002 


-0.06 


7.1219E-004 


0.94 


Grid 


llr - T h II 

1 1 '-cy xy\\ oo 


Order 


\\Txv -T* 111 
ll y xy 1 1 J- 


Order 


16 2 


1.6864E-001 




8.7332E-002 




32 2 


1.0266E-001 


0.72 


3.9547E-002 


1.14 


64 2 


2.1950E-001 


-1.10 


1.8694E-002 


1.08 


128 2 


2.0094E-001 


0.13 


6.9612E-003 


1.43 


256 2 


2.6596E-001 


-0.40 


4.0559E-003 


0.78 


512 2 


2.0961E-001 


0.34 


1.6447E-003 


1.30 


1024 2 


2.5255E-001 


-0.27 


8.8242E-004 


0.90 


Grid 


\\u- U h \\oo 


Order 


||u-it ft ||i 


Order 


16 2 


3.4171E-001 




1.9727E-001 




32 2 


7.4246E-002 


2.20 


4.3033E-002 


2.20 


64 2 


2.6593E-002 


1.48 


1.2321E-002 


1.80 


128 2 


9.2292E-003 


1.53 


3.3497E-003 


1.88 


256 2 


6.7182E-003 


0.46 


1.5327E-003 


1.13 


512 2 


3.0843E-003 


1.12 


5.9129E-004 


1.37 


1024 2 


1.7877E-003 


0.79 


2.7579E-004 


1.10 



zero traction condition (|3.5|) enforces a relationship between the surface pressure and 
the viscous stress resulting from this velocity field. To satisfy this condition, we use 
the following expression for pressure: 

p = 0l l^ r 2 pt sin(2fl) sin(%/3 In r) (15 - 30r + 16r 2 ) (9.2) 
81 

The pressure in this expression will be non-zero at the interface; any method to solve 
this problem will need to correctly handle the coupling between pressure and viscous 
stresses. From this information, the expressions for the input velocity and final stresses 
can be derived using equations (|4.1[) - (|4.3p . We used a computer algebra system for 
this purpose. The convergence results are shown in Table 19.21 and Figure 19.11 

9.2. Stokes Flow with Solid Wall Boundaries (2D). Our Stokes solid 
boundary test case is an annulus centred at the origin with inner radius r = 0.5, 



FINITE DIFFERENCE STOKES FLOW ON IRREGULAR DOMAINS 



Pressure, L 



Pressure, LI 




o 



-1.2- 
-1.4 



1 



1.5 2 2.5 3 3.5 

loglO(N), for grid size NxN 

(a) 

Stress x , L°° 




1.5 2 2.5 3 

log 1 0(N), for grid size NxN 



-0.5 

-1 
-1.5 

-2- 
-2.5 

-3 
-3.5 



1 



1.5 2 2.5 3 

loglO(N), for grid size NxN 

(b) 

Stress x , LI 




1.5 2 2.5 3 3.5 

loglO(N), for grid size NxN 



(c) 



(d) 



Stress x , L« 



Stress x , LI 

xy' 




1.5 2 2.5 3 

log 1 0(N), for grid size NxN 

(g) 



-1 
-1.5 

-2 
-2.5- 

-3 : 
-3.5 



-0.5 n 

-1 
-1.5 

-2 
-2.5 

-3 
-3.5 

-4 



1 




1.5 2 2.5 3 3.5 

loglO(N), for grid size NxN 



(f) 

Velocity, LI 




1.5 2 2.5 3 3.5 

loglO(N), for grid size NxN 



Fig. 9.1. Convergence graphs for the Stokes problem with free surface boundaries (2D). 
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Table 9.3 

Convergence of Stokes with solid walls (2D) 



Grid 


Hp - P /l llno 

\\f r || oo 


Order 


\\r y 111 


Order 


16 2 


3.4118E+000 




2.1013E+000 




32 2 


3.6933E+000 


-0.11 


9.0009E-001 


1.22 


64 2 


3.0338E+000 


0.28 


4.5162E-001 


0.99 


128 2 


2.5942E+000 


0.23 


1.5373E-001 


1.55 


256 2 


2.4222E+000 


0.10 


8.5695E-002 


0.84 


512 2 


2.3573E+000 


0.04 


3.9879E-002 


1.10 


1024 2 


2.3381E+000 


0.01 


2.0266E-002 


0.98 


Grid 


Ik - T h II 
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Order 
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Order 
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Order 
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0.97 



outer radius r = 1, density p = 1, and viscosity [i = 0.1, computed over a timestep 
At = 1. Inner and outer boundaries are static solids. We will again use a streamfunc- 
tion ip to dictate our velocity field and ensure it is divergence free: 

= 256r 4 - 768r 3 + 832r 2 - 384r + 64 (9.3) 

For pressure, we use: 

p = r 2 cos(6>) sin(0) (9.4) 

Equations (14. f)4.3|) can be used to derive the input velocities and final stresses. The 
convergence results are shown in Table 19.31 and Figure 19.21 

9.3. Stokes Flow with Both Solid Wall and Free Surface Boundaries 
(2D). To test a scenario featuring both solid and free surface boundaries, we solve 



FINITE DIFFERENCE STOKES FLOW ON IRREGULAR DOMAINS 



Pressure, L>° 



Pressure, LI 



0.7 
0.65 
0.6 
O 0.55 



0.4 
0.35 



0.2 
0.15 



2 o.i 



O) 0.05- 
O 

0- 
-0.05 



■ Error 
- 0th order guide 




1.5 2 2.5 3 

loglO(N), for grid size NxN 

(a) 

Stress x , L~ 




1.5 2 2.5 3 

loglO(N), for grid size NxN 



o- 

-0.5- 



-1.5- 
-2 




0.5- 
0- 
-0.5- 

-1 - 
-1.5 : 
-2 



1 



1.5 2 2.5 3 

log10(N), for grid size NxN 

(b) 

Stress x , LI 




1.5 2 2.5 3 3.5 

loglO(N), for grid size NxN 



(c) 



(d) 



0.5 
0.4 



2 0.3 



a, 0.2 
o 




Stress x LI 



O)-0.5 
O 



1.5 2 2.5 3 3.5 

log 1 0(N), for grid size NxN 

(e) 

Velocity, L°° 





-•-Error 

1st order guide 








1 1.5 2 


2.5 3 3 


loglO(N), for grid size NxN 


(g) 





-1.5 
-2 



1 

0.5- 


-0.5 
-1 

-1.5 
-2- 

-2.5- 





-•-Error 




1st order guide 















1.5 2 2.5 3 3.5 

log10(N), for grid size NxN 

(f) 

Velocity, LI 





Error 




1st order guide 











1.5 2 2.5 3 

log10(N), for grid size NxN 



Fig. 9.2. Convergence graphs for the Stokes problem with solid wall boundaries (2D). 
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for fluid motion in an annulus centred at the origin with inner radius r = 0.1, outer 
radius r = 0.75, density p = 1, and viscosity fi — 0.1, over a timestep At — 1. The 
outer boundary is a free surface, and the inner boundary is a static solid. We will 
again use a streamfunction ip to dictate our velocity held and ensure it is divergence 
free: 

ip = o 8 _ 1 o n ° o cos(20)r 4 cos(%/3lnr)(169 - 390r + 240r 2 ) (9.5) 
371293 

For pressure, we use: 

320000 



P = 



371293 



VSfir 2 sin(20) sin(%/31nr)(169 - 390r + 240r 2 ) (9.6) 

Equations (14. ()4.3|) can be used to derive the input velocities and final stresses. The 
convergence results are shown in Table 19.41 and Figure [ 



9.4. Stokes Flow with Prescribed Velocity Solid Boundaries (2D). To 

test solid boundaries with prescribed (non-zero) boundary velocities, we solve for fluid 
motion in an annulus centred at the origin, with inner radius r — 0.5, outer radius 
r = 1, density p = 1, and viscosity /i = 0.1, over a timestep At = 1. The outer 
boundary is static, while the inner boundary rotates rigidly with a clockwise angular 
velocity lo = 2. For the final velocity we use the streamfunction ip: 



For pressure, we use: 



tb = r 4 - 3r 3 + ^r 2 + ]-r + \ (9.7) 

4 2 4 y ' 



p = r 2 cos(0) sin(0) (9.8) 



As in the preceding examples, stresses and input velocities can be computed from 
equations (|4.1l) - (|4.3p . Convergence results are shown in Table l9~5l and Figure! 



9.5. Three Dimensional Stokes Flow. To test our method in three dimen- 
sions where analytical solutions are substantially more difficult to derive, we created 
a numerical solution at resolution 155 3 , and tested convergence towards this solution. 
The test case consists of a sphere of liquid centred at the origin with a free surface 
at r = 1, containing a nested static solid sphere of radius r = 0.25. Density was set 
to p = 1 and viscosity to p, = 0.1. We compute a time step of length At = 1 starting 
from an input velocity field: 



>inp 



ut = (0.5 + x sin(yz), cos(0.25y) + 0.5xyz, x + Q.5yz 2 + \Jl + y 2 ) (9.9) 



Convergence results are shown in Table 19.61 Interestingly, this numerical experiment 
suggests first order convergence in L°° even for stress variables, an improvement over 
the results in 2D. We cannot yet explain why the move to 3D would bring greater 
accuracy, but several similar tests showed the same pattern. 



9.6. Discussion. Table 19771 summarizes the approximate orders of convergence 
suggested by our experiments in two dimensions. As noted earlier, we achieve essen- 
tially first order convergence, with the exception of pressure and stress in L°° . We 
now proceed to make some general comments. 

A rigourous convergence theory to describe the method is beyond the scope of the 
current work. However, in the absence of boundaries, our methods are equivalent to 
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Fig. 9.3. Convergence graphs for the Stokes problem with both a free surface and a solid wall 
boundary (2D). 



FINITE DIFFERENCE STOKES FLOW ON IRREGULAR DOMAINS 



21 




Fig. 9.4. Convergence graphs for the Stokes problem with a moving solid wall boundary (2D). 
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Table 9.4 

Convergence of Stokes with solid walls and a free surface (2D) 



Grid 


Hp - p'MIoo 

WJr' I II OO 


Order 


Ilp-P ft lh 
\\y v 1 1 -L 


Order 


16 2 


2.2404E-001 




2.1048E-001 




32 2 


2.1216E-001 


0.08 


6.5787E-002 


1.68 


64 2 


1.6222E-001 


0.39 


2.2932E-002 


1.52 


128 2 


1.2874E-001 


0.33 


7.8514E-003 


1.55 


256 2 


1.7168E-001 


-0.42 


5.1902E-003 


0.60 


512 2 


1.7212E-001 


-0.00 


2.1723E-003 


1.26 


1024 2 


1.6865E-001 


0.03 


1.1268E-003 


0.95 


Grid 


Ht - T h II 

1 1 xx XX 1 1 °° 


Order 


| T rT — 7"'* | 1 


Order 


16 2 


4.3812E-001 




2.0762E-001 




32 2 


1.8386E-001 


1.25 


5.7557E-002 


1.85 


64 2 


9.5150E-002 


0.95 


1.6510E-002 


1.80 


128 2 


4.9719E-002 


0.94 


5.9357E-003 


1.48 


256 2 


6.7733E-002 


-0.45 


3.3888E-003 


0.81 


512 2 


6.9916E-002 


-0.05 


1.3958E-003 


1.28 


1024 2 


7.3807E-002 


-0.08 


7.2576E-004 


0.94 


Grid 


Ht - T h II 

1 1 xy xy 1 1 oo 


Order 


Iheu-T* 111 

ll J- y xy 1 1 J- 


Order 


16 2 


2.1174E-001 




1.0047E-001 




32 2 


1.2921E-001 


0.71 


4.8916E-002 


1.04 


64 2 


2.6177E-001 


-1.02 


2.0869E-002 


1.23 


128 2 


2.1197E-001 


0.30 


7.1786E-003 


1.54 


256 2 


2.7070E-001 


-0.35 


3.9133E-003 


0.88 


512 2 


2.1190E-001 


0.35 


1.5764E-003 


1.31 


1024 2 


2.5361E-001 


-0.26 


8.5715E-004 


0.88 


Grid 


\\u - U h \\oo 


Order 




Order 


16 2 


4.2290E-001 




2.3654E-001 




32 2 


8.6887E-002 


2.28 


4.8516E-002 


2.29 


64 2 


2.9712E-002 


1.55 


1.3486E-002 


1.85 


128 2 


1.0943E-002 


1.44 


3.6100E-003 


1.90 


256 2 


5.3552E-003 


1.03 


1.3797E-003 


1.39 


512 2 


2.8869E-003 


0.89 


5.2047E-004 


1.41 


1024 2 


1.6769E-003 


0.78 


2.4870E-004 


1.07 



a straightforward discretization of the original PDE form with second order centred 
finite differences on a staggered grid, and can therefore expect to achieve uniform 
second order convergence. Near boundaries this clearly does not hold, and the effective 
quadrature is likely only first order due to the use of piecewise constant approximations 
of integral terms. This is in line with our results in L°° for velocity. Solution gradients 
can generally be expected to be one order less accurate than the solution itself, and 
this is also evident in the O(l) errors in the L°° results for fluid pressure and stress. 
Interestingly, Chen et al. recently showed that the immersed boundary method for 
Stokes flow likewise exhibits 0(1) errors in pressure in L°° [5]. 

The virtual node method of Bcdrossian et al. [4] for the Poisson equation uses a 
variational formulation similar to ours, but makes use of piecewise bilinear Cartesian 
elements near the boundary to estimate the relevant integrals, at the cost of denser 
stencils for boundary cells. Their results indicate second order convergence which is 
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Table 9.5 

Convergence of Stokes with prescribed velocity solid walls (2D) 



Grid 


Hp- p^IIoo 

\\F f || oo 


Order 


\\r y 111 


Order 


16 2 


4.7210E-002 




2.2007E-002 




32 2 


2.7144E-002 


0.80 


7.9999E-003 


1.46 


64 2 


5.7009E-002 


-1.07 


5.0198E-003 


0.67 


128 2 


4.5540E-002 


0.32 


2.2860E-003 


1.13 


256 2 


5.5697E-002 


-0.29 


1.3316E-003 


0.78 


512 2 


6.0562E-002 


-0.12 


6.4669E-004 


1.04 


1024 2 


6.1643E-002 


-0.03 


3.3755E-004 


0.94 


Grid 


Ik . - T h II 

1 1 X X i ' 1 1 OC 


Order 


| r TT — riLlh 

II ■L-L XX II 1 


Order 


16 2 


2.6730E-002 




1.2788E-002 




32 2 


1.1167E-002 


1.26 


6.2639E-003 


1.03 


64 2 


2.2172E-002 


-0.99 


4.4126E-003 


0.51 


128 2 


2.0462E-002 


0.12 


1.8498E-003 


1.25 


256 2 


3.0896E-002 


-0.59 


1.0722E-003 


0.79 


512 2 


2.7147E-002 


0.19 


5.2769E-004 


1.02 


1024 2 


2.8528E-002 


-0.07 


2.7240E-004 


0.95 


Grid 


llr - T h II 

1 1 '£y xy 1 1 oo 


Order 


\\t xv -t* 111 

ll y xy 1 1 J- 


Order 


16 2 


1.8177E-002 




1.6688E-002 




32 2 


1.9415E-002 


-0.10 


8.3726E-003 


1.00 


64 2 


2.4280E-002 


-0.32 


4.4206E-003 


0.92 


128 2 


2.6202E-002 


-0.11 


2.0984E-003 


1.07 


256 2 


3.2232E-002 


-0.30 


1.1682E-003 


0.85 


512 2 


3.3983E-002 


-0.08 


5.9705E-004 


0.97 


1024 2 


3.3056E-002 


0.04 


3.0660E-004 


0.96 


Grid 


\\u- U h \\oo 


Order 


||u-it ft ||i 


Order 


16 2 


6.6980E-002 




2.9655E-002 




32 2 


4.2451E-002 


0.66 


9.2806E-003 


1.68 


64 2 


1.9403E-002 


1.13 


2.5526E-003 


1.86 


128 2 


8.5696E-003 


1.18 


8.4643E-004 


1.59 


256 2 


3.7606E-003 


1.19 


3.5600E-004 


1.25 


512 2 


2.4402E-003 


0.62 


1.4966E-004 


1.25 


1024 2 


1.4085E-003 


0.79 


7.1476E-005 


1.07 



consistent with the fact that our use of piecewise constant estimates yields first order 
convergence. This also suggests that applying bilinear elements near boundaries may 
be effective in raising the convergence order of our method for Stokes flow, while 
maintaining the benefits of sparsity and positive-definiteness. 

Ng et al. 33 j pointed out that in the final discretized form, their method for the 
Poisson equation with Neumann solid boundaries is identical to that of Batty et al. 
[2] if the face volume weights suggested by the variational perspective are replaced 
with face area (finite volume) weights. The latter choice leads to an increase in L°° 
accuracy for pressure from first to second order and velocity from zeroth to first order. 
The related ghost fluid method for the Poisson equation with Dirichlet boundaries [20] 
likewise exhibits second order in pressure and first order in velocity, with a different 
choice of weights. While this hints that alternative diagonal weighting matrices W 
might raise the order of accuracy of the current method, we have found that directly 
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Table 9.6 

Convergence of 3D Stokes with solid walls and a free surface 



Grid 


\\p-p h \\oo Order 


p — p h i Order 


12 2 
24 2 
48 2 


1.3912E-002 
7.5008E-003 0.89 
4.9163E-003 0.61 


1.3312E-002 
6.3879E-003 1.06 
3.5866E-003 0.83 


Grid 


\\t X x ~ t^Woc Order 


\\l~xx — TxxW^ 0rdcr 


12 2 

24 2 
48 2 


1.7478E-002 
1.0965E-002 0.67 
5.8661E-003 0.90 


2.0401E-002 
1.1213E-002 0.86 
5.6263E-003 0.99 


Grid 


\\r yy - t^Joo Order 


\\ T VV- T yy\\l 0rdcr 


12 2 

24 2 
48 2 


2.3604E-002 
1.3516E-002 0.80 
6.8033E-003 0.99 


2.9434E-002 
1.5741E-002 0.90 
7.7348E-003 1.03 


Grid 


\\ T xy - Tx V \\oo Order 


hxy-r^h Order 


12 2 

24 2 
48 2 


6.4309E-003 
5.1895E-003 0.31 
3.1232E-003 0.73 


7.9982E-003 
4.4602E-003 0.84 
2.0025E-003 1.16 


Grid 


Ikxz-T^Hoo Order 


Wr^-r^h Order 


12 2 

24 2 
48 2 


5.5459E-002 
3.0026E-002 0.89 
1.3128E-002 1.19 


1.3162E-001 
7.6418E-002 0.78 
3.6082E-002 1.08 


Grid 


\\t V z ~ rfzWoo Order 


Wryz-r^h Order 


12 2 
24 2 
48 2 


1.6390E-002 
1.0325E-002 0.67 
5.2896E-003 0.96 


2.8919E-002 
1.7670E-002 0.71 
8.7857E-003 1.01 


Grid 


\\u — u h oo Order 


\\u — u h \\i Order 


12 2 

24 2 
48 2 


1.9163E-001 
1.2664E-001 0.60 
7.1250E-002 0.83 


2.5741E-001 
1.4041E-001 0.87 
6.3736E-002 1.14 


Grid 


\\v — v h oo Order 


\\v — v h \\i Order 


12 2 

24 2 
48 2 


6.2893E-002 
4.7289E-002 0.41 
2.9547E-002 0.68 


5.4328E-002 
3.6334E-002 0.58 
2.0160E-002 0.85 


Grid 


\\w — w h oo Order 


\\w — w h i Order 


12 2 
24 2 
48 2 


1.8633E-001 
1.2850E-001 0.54 
7.4279E-002 0.79 


2.7837E-001 
1.5601E-001 0.84 
7.1820E-002 1.12 



introducing ghost fluid or finite volume weights in this setting breaks the symmetry 
of the linear system. Nonetheless, a deeper exploration of the connections between 
our method and ghost fluid/immersed interface methods, finite volume methods, and 
finite element methods might provide the key to an improved weighting scheme. 

10. Application to Viscous Jet Buckling. One particularly fascinating phe- 
nomenon exhibited by highly viscous liquids is jet buckling. When a falling liquid 
column of sufficient viscosity impacts a solid surface, it will fold or coil over on itself 
rather than spreading out smoothly. Relatively few researchers have looked at sim- 
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Table 9.7 
Convergence Behaviour in 2D 



Variable 


L°° 


L 1 


Pressure 





1 


Deviatoric Stress 





1 


Velocity 


1 


1 



ulating Newtonian viscous buckling, despite its prevalence in many common liquids 
such as honey. To the best of our knowledge, the GENSMAC code of Tome, McKee 
and co-authors is the only prior finite difference scheme to do so in three dimensions 
[2HJ [501 EHl [35] . However, as noted earlier this approach requires a case-by-case anal- 
ysis of discrete surface orientations and its implicit formulation entails solving a large 
non-symmetric linear system. Jet buckling has also been addressed in a finite element 
setting [5] and with an SPH approach [4"T] . 

With this problem in mind, we incorporate our Stokes solver into a simple two- 
stage fractional step Navier-Stokes routine (similar to that presented in section 2 of 
the paper by Ng et al. [33]). First, starting with a velocity field it 71 at time t n , we 
compute advection and body forces to produce an intermediate velocity field u*: 



At 



(VtT) = F (10.1) 



We account for advection terms with a first order semi-Lagrangian scheme using 
bilinear interpolation of velocities. We then simply add any external body forces F 
(gravity in our examples). From this intermediate velocity, we then simultaneously 
incorporate viscous forces and project the velocity field to be divergence free using 
our Stokes solver, to arrive at time t n+1 = t n + At: 

p [ ' = V • r" +1 - Vp n+1 (10.2) 

V-m = (10.3) 
r" +1 = M (V?T +1 + (W" +1 ) T ) (10.4) 

with the appropriate free surface and solid boundary conditions applied. Tracking 
of the liquid surface position is performed using a basic semi-Lagrangian level set 
method (eg. [16]). 

10.1. Two Dimensional Jet Buckling. Figure [TO . 1 1 presents the results of a 
two-dimensional simulation of planar viscous jet buckling. The simulation domain is a 
circle of radius 0.4[to] centred at (0.5, 0.5). A horizontal ceiling is placed at y = 0.8[m], 
featuring a liquid jet inflow centered at x = 0.5 [to] with a fixed vertical velocity of 
U = — 0.5[to/s] and a width of D — 0.06[m]. This configuration yields a drop height 
of H = 0.7[m]. The density of the liquid is p — l[kg/m 3 ] and the dynamic viscosity 
of the liquid is /i = 0.075[Pa ■ s]. Gravity is set at — 9.81[m/s 2 ]. The simulation grid 
used a resolution of 150 x 150 cells. 

Following Tome and McKee [50], this yields a Reynolds number of Re = = 
0.4 and an aspect ratio for the liquid jet of H/D = 0.7/0.06 = 11.667. This falls 
within the regime in which planar buckling is expected to occur (Re < 0.5, H/D > 10) 
according to Cruikshank and MunsonflTl |10j; as shown our method reproduces the 
buckling phenomenon. 




Fig. 10.1. A two-dimensional example of viscous jet buckling performed using our simple 
Navier-Stokes routine. The first image occurs 0.5 seconds into the simulation, and subsequent 
frames occur at 0.2 second intervals. 
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10.2. Three Dimensional Jet Buckling. Figures ri0.2l and [10.3l present the re- 
sults of a three-dimensional simulation of cylindrical viscous jet buckling (ie. coiling). 
The simulation domain is a sphere of radius 0.4 [to] centred at (0.5, 0.5, 0.5). A circu- 
lar inlet is centred at (0.5,0.8,0.5), with a fixed vertical velocity of U = — 0.5[to/s] 
and a diameter D — 0.08[to]. This configuration yields a drop height of H = 0.7[m]. 
The density of the liquid is p = l[kg/m 3 ] and the dynamic viscosity of the liquid is 
fj, = 0.3[Pa ■ s]. Gravity is set at — 9.81[to/s 2 ]. The simulation grid used a resolution 
of 80 x 80 x 80 cells. 

This yields a Reynolds number of Re = k, 0.133 and an aspect ratio for 

the liquid jet of H/D = 0.7/0.08 = 8.75. This falls within the guidelines for when 
axisymmetric buckling typically arises (Re < 1.2, H/D > 7) according to Cruikshank 
and Munson[IT1 IIP], and the result does indeed exhibit substantial buckling. 

11. Conclusions and Future Work. We have shown that a Cartesian grid 
finite difference method derived from a variational principle can correctly capture dif- 
ficult irregular boundary conditions in Stokes flow problems, while providing stability 
for large time steps and yielding a sparse, symmetric positive-definite linear system. 
To do so, we have unified and extended recent work on embedded boundary meth- 
ods for pressure projection and viscosity. In practice the method's implementation 
is remarkably simple, yet it correctly captures the challenging free surface bound- 
ary condition that enables simulation of jet buckling phenomena. This work raises a 
number of questions and directions for future work. 

In our numerical study, we observed improved L°° convergence in 3D compared 
to 2D, and plan to investigate if this truly holds — perhaps beginning by deriving a 
full-fledged analytic test case as we have done in 2D. The 2D convergence test cases we 
presented also considered only scenarios where the two different boundary conditions 
(solid and free surface) do not meet. We suspect that this is an inherently more 
difficult problem to address, giving rise to issues analogous to those which occur in 
the presence of sharp boundary features, and our preliminary experiments support 
this conjecture. Nevertheless, such configurations occur frequently in the buckling 
examples we have included, illustrating that the method remains stable and provides 
qualitatively reasonable results. 

While the viscous jet buckling example provides a practical validation of our 
method's boundary condition enforcement, the current underlying Navier-Stokes sim- 
ulator is fairly basic. A thorough study of this phenomenon would likely need to 
consider improved advection and time-splitting methods in place of the first order 
approaches applied here. 

In terms of handling related phenomena, we have noted that our work is closely 
related to that of Batty et al. who considered considered the simpler Poisson problem 
for incompressibility with rigid bodies [2]. An obvious extension of the current work 
would therefore be to consider Stokes flow coupled with fully dynamic deformable 
structures. Two-phase flow would also be a useful direction to pursue, as would non- 
Newtonian fluid models. Finally, it would be interesting to consider whether exploiting 
variational principles in this manner might be useful for handling irregular boundaries 
in other problems that are commonly discretized on staggered grids. Some potential 
examples include vorticity-based formulations of fluid flow, diffusion problems, or 
porous flow. 

Appendix A. Notation. 

The following symbols and letters are used throughout the paper, with units given 
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(g) (h) (i) 

Fig. 10.2. A three-dimensional example of viscous jet buckling performed using our simple 
Navier-Stokes routine. The first image occurs 0.6 seconds into the simulation, and subsequent 
frames occur at 0.3 second intervals. Additional images are shown in Figure D 0. 3\ 



for dimensionful quantities: 

fi : coefficient of dynamic viscosity, Pa ■ s 
p : density coefficient, kg/m 3 
flp : fluid domain 

51s : solid domain (the complement of CIf) 
51l : liquid domain 

Q a ■ air domain (the complement of Ql) 

u : velocity vector, m/s 

p : pressure, Pa 

t : deviatoric stress tensor, Pa 
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(g) (h) (i) 



Fig. 10.3. Additional images of viscous jet buckling, continued from Figure D 0. 2\ 

At : time step, s 

T : traction vector, Pa 

D : discrete deformation rate matrix 

G : discrete gradient matrix 

M : diagonal matrix of viscosity coefficients, per velocity sample 

P : diagonal matrix of density coefficients, per pressure/stress sample 

Wf '■ diagonal matrix of fluid (non-solid) fraction weights 

Ws ■ diagonal matrix of solid fraction weights, complementing Wl 

Wl : diagonal matrix of liquid (non-air) fraction weights 

Wa ■ diagonal matrix of air fraction weights, complementing Wl 

W u , W p , W T : superscripts on weight matrices indicate the associated sample 
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position, ie. velocity u (cell faces), pressure p (cell centres), stresses r (cell 
centres and edges). 
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